//该文件中，包含了普通的MPI优化，优化后的MPI优化，循环划分和块划分的MPI优化
//非阻塞通信的MPI优化
#include<iostream>
#include <arm_neon.h>
#include<stdio.h>
#include<typeinfo>
#include<arm_neon.h>
#include<stdlib.h>
#include<cmath>
#include<mpi.h>
#include<string.h>
#include <omp.h>
using namespace std;
int N = 2000;
#define NUM_THREADS 7
float** A = NULL;

void A_init() 
{
    A = new float *[N];
    for(int i=0; i<N; i++)
    {
        A[i] = new float[N];
    }
    for(int i=0;i<N;i++)
    {
        for(int j=0; j<i; j++)
        {
            A[i][j]=0;
        }
        A[i][i]=1.0;
        for(int j=i+1;j<N;j++)
        {
            A[i][j]=static_cast <float> (rand()) / (static_cast <float> (RAND_MAX/1000.0));
        }
    }
    for(int k=0;k<N;k++)
        for(int i=k+1;i<N;i++)
            for(int j=0;j<N;j++)
            {
                A[i][j]+=A[k][j];
            }
}

void A_initAsEmpty() 
{
    A = new float *[N];
    for(int i=0; i<N; i++)
    {
        A[i] = new float[N];
    }
    for(int i=0; i<N; i++)
    {
        for(int j=0; j<N; j++)
        {
            A[i][j] = 0;
        }
    }
}

void deleteA() 
{
    for (int i = 0; i < N; i++) 
    {
        delete[] A[i];
    }
    delete A;
}

void print(float** a) 
{
    for (int i = 0; i < N; i++) 
    {
        for (int j = 0; j < N; j++) 
        {
            cout << a[i][j] << " ";
        }
        cout << endl;
    }
}

void xiaoyuan_cx() 
{
    double start_time = 0;
    double end_time = 0;
    start_time = MPI_Wtime();
    A_init();
    for (int k = 0; k < N; k++) 
    {
        for (int j = k + 1; j < N; j++) 
        {
            A[k][j] = A[k][j] / A[k][k];
        }
        A[k][k] = 1.0;
        for (int i = k + 1; i < N; i++) {
            for (int j = k + 1; j < N; j++) {
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            }
            A[i][k] = 0;
        }
    }
    end_time = MPI_Wtime();
    printf("串行算法耗时:%.4lf ms\n", 1000 * (end_time - start_time));
}

void mpi_talk(int argc, char* argv[])
{
    double start_time = 0;
    double end_time = 0;
    int total = 0;
    int rank = 0;
    int i = 0;
    int j = 0;
    int k = 0;
    MPI_Status status;
    MPI_Comm_size(MPI_COMM_WORLD, &total);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    int begin = N / total * rank;
    int end = (rank == total - 1) ? N : N / total * (rank + 1);
    if (rank == 0) 
    {
        A_init();
        //print(A);
        //cout << endl << endl;
        start_time = MPI_Wtime();
        for (j = 1; j < total; j++) 
        {
            int b = j * (N / total), e = (j == total - 1) ? N : (j + 1) * (N / total);
            for (i = b; i < e; i++) 
            {
                MPI_Send(&A[i][0], N, MPI_FLOAT, j, 1, MPI_COMM_WORLD);//1是初始矩阵信息，向每个进程发送数据
            }
        }
    }
    else 
    {
        A_initAsEmpty();
        for (i = begin; i < end; i++) 
        {
            MPI_Recv(&A[i][0], N, MPI_FLOAT, 0, 1, MPI_COMM_WORLD, &status);
        }

    }

    MPI_Barrier(MPI_COMM_WORLD);  //此时每个进程都拿到了数据
    //start_time = MPI_Wtime();
    for (k = 0; k < N; k++) //负责第 k 行的进程进行除法运算，将除法结果一对多广播给编号更大的进程
    {
        if ((begin <= k && k < end)) 
        {
            for (j = k + 1; j < N; j++) 
            {
                //A[k][j] = A[k][j] / A[k][k];
            }
            //A[k][k] = 1.0;
            for (j = 0; j < total; j++) 
            {
                if(j!=rank)
                    MPI_Send(&A[k][0], N, MPI_FLOAT, j, 0, MPI_COMM_WORLD);//0号消息表示除法完毕
            }
        }
        else
        {
            int src;
            if (k < N / total * total)//在可均分的任务量内
                src = k / (N / total);
            else
                src = total-1;
            MPI_Recv(&A[k][0], N, MPI_FLOAT, src, 0, MPI_COMM_WORLD, &status);
        }
        for (i = max(begin, k + 1); i < end; i++) 
        {
            for (j = k + 1; j < N; j++) 
            {
                //A[i][j] = A[i][j] - A[i][k] * A[k][j];
            }
            //A[i][k] = 0;
        }
    }
    MPI_Barrier(MPI_COMM_WORLD);	//各进程同步
    if (rank == 0) 
    {//0号进程中存有最终结果
        end_time = MPI_Wtime();
        printf("MPI(纯通信)耗时:%.4lf ms\n", 1000 * (end_time - start_time));
        //print(A);
    }
}

void xiaoyuan_mpi(int argc, char* argv[]) //块划分 
{
    double start_time = 0;
    double end_time = 0;
    int total = 0;
    int rank = 0;
    int i = 0;
    int j = 0;
    int k = 0;
    MPI_Status status;
    MPI_Comm_size(MPI_COMM_WORLD, &total);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    int begin = N / total * rank;
    int end = (rank == total - 1) ? N : N / total * (rank + 1);
    if (rank == 0) 
    {
        A_init();
        start_time = MPI_Wtime();
        for (j = 1; j < total; j++) 
        {
            int b = j * (N / total), e = (j == total - 1) ? N : (j + 1) * (N / total);
            for (i = b; i < e; i++) 
            {
                MPI_Send(&A[i][0], N, MPI_FLOAT, j, 1, MPI_COMM_WORLD);
            }
        }
    }
    else 
    {
        A_initAsEmpty();
        for (i = begin; i < end; i++) 
        {
            MPI_Recv(&A[i][0], N, MPI_FLOAT, 0, 1, MPI_COMM_WORLD, &status);
        }

    }

    MPI_Barrier(MPI_COMM_WORLD);
    for (k = 0; k < N; k++)
    {
        if ((begin <= k && k < end)) 
        {
            for (j = k + 1; j < N; j++) 
            {
                A[k][j] = A[k][j] / A[k][k];
            }
            A[k][k] = 1.0;
            for (j = 0; j < total; j++) 
            {
                if(j!=rank)
                    MPI_Send(&A[k][0], N, MPI_FLOAT, j, 0, MPI_COMM_WORLD);
            }
        }
        else
        {
            int src;
            if (k < N / total * total)
                src = k / (N / total);
            else
                src = total-1;
            MPI_Recv(&A[k][0], N, MPI_FLOAT, src, 0, MPI_COMM_WORLD, &status);
        }
        for (i = max(begin, k + 1); i < end; i++) 
        {
            for (j = k + 1; j < N; j++) 
            {
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            }
            A[i][k] = 0;
        }
    }
    MPI_Barrier(MPI_COMM_WORLD);
    if (rank == 0) 
    {
        end_time = MPI_Wtime();
        printf("MPI块划分耗时:%.4lf ms\n", 1000 * (end_time - start_time));
    }
}

void xiaoyuan_mpi_simd(int argc, char* argv[]) //块划分+simd
{
    double start_time = 0;
    double end_time = 0;
    int total = 0;
    int rank = 0;
    int i = 0;
    int j = 0;
    int k = 0;
    MPI_Status status;
    MPI_Comm_size(MPI_COMM_WORLD, &total);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    //cout << rank << " of " << total << " created" << endl;
    int begin = N / total * rank;
    int end = (rank == total - 1) ? N : N / total * (rank + 1);
    //cout << "rank " << rank << " from " << begin << " to " << end << endl;
    if (rank == 0) 
    {  //0号进程初始化矩阵
        A_init();
        //cout << "initialize success:" << endl;
        //print(A);
        //cout << endl << endl;
        start_time = MPI_Wtime();
        for (j = 1; j < total; j++) 
        {
            int b = j * (N / total), e = (j == total - 1) ? N : (j + 1) * (N / total);
            for (i = b; i < e; i++) 
            {
                MPI_Send(&A[i][0], N, MPI_FLOAT, j, 1, MPI_COMM_WORLD);//1是初始矩阵信息，向每个进程发送数据
            }
        }
    }
    else 
    {
        A_initAsEmpty();
        for (i = begin; i < end; i++) 
        {
            MPI_Recv(&A[i][0], N, MPI_FLOAT, 0, 1, MPI_COMM_WORLD, &status);
        }

    }

    MPI_Barrier(MPI_COMM_WORLD);  //此时每个进程都拿到了数据
    //start_time = MPI_Wtime();
    for (k = 0; k < N; k++) //负责第 k 行的进程进行除法运算，将除法结果一对多广播给编号更大的进程
    {
        if ((begin <= k && k < end)) 
        {
            for (j = k + 1; j < N; j++) 
            {
                A[k][j] = A[k][j] / A[k][k];
            }
            A[k][k] = 1.0;
            for (j = 0; j < total; j++)
            {
                if(j!=rank)
                    MPI_Send(&A[k][0], N, MPI_FLOAT, j, 0, MPI_COMM_WORLD);//0号消息表示除法完毕
            }
        }
        else
        {
            int src;
            if (k < N / total * total)//在可均分的任务量内
                src = k / (N / total);
            else
                src = total-1;
            MPI_Recv(&A[k][0], N, MPI_FLOAT, src, 0, MPI_COMM_WORLD, &status);
        }
        for (i = max(begin, k + 1); i < end; i++) 
        {
            float32x4_t vaik;
            vaik = vdupq_n_f32(A[i][k]);
            int j;
            for (j = k + 1; j + 4 <= N; j += 4) 
            {
                float32x4_t vakj, vaij, vx;
                vakj = vld1q_f32(&A[k][j]);
                vaij = vld1q_f32(&A[i][j]);
                vx = vmulq_f32(vakj, vaik);
                vaij = vsubq_f32(vaij, vx);
                vst1q_f32(&A[i][j], vaij);
            }
            for (; j < N; j++) 
            {
                A[i][j] = A[i][j] - A[k][j] * A[i][k];
            }
            A[i][k] = 0;
        }
    }
    MPI_Barrier(MPI_COMM_WORLD);	//各进程同步
    if (rank == 0) 
    {//0号进程中存有最终结果
        end_time = MPI_Wtime();
        printf("MPI块划分+simd耗时:%.4lf ms\n", 1000 * (end_time - start_time));
        //print(A);
    }
}

void xiaoyuan_mpi_openmp(int argc, char* argv[]) //块划分+openmp 
{
    double start_time = 0;
    double end_time = 0;
    int total = 0;
    int rank = 0;
    int i = 0;
    int j = 0;
    int k = 0;
    MPI_Status status; //*
    MPI_Comm_size(MPI_COMM_WORLD, &total);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    //cout << rank << " of " << total << " created" << endl;
    int begin = N / total * rank;
    int end = (rank == total - 1) ? N : N / total * (rank + 1);
    //cout << "rank " << rank << " from " << begin << " to " << end << endl;
    if (rank == 0) 
    {  //0号进程初始化矩阵
        A_init();
        //cout << "initialize success:" << endl;
        //print(A);
        //cout << endl << endl;
        start_time = MPI_Wtime();
        for (j = 1; j < total; j++) 
        {
            int b = j * (N / total), e = (j == total - 1) ? N : (j + 1) * (N / total);
            for (i = b; i < e; i++) 
            {
                MPI_Send(&A[i][0], N, MPI_FLOAT, j, 1, MPI_COMM_WORLD);//1是初始矩阵信息，向每个进程发送数据
            }
        }
    }
    else 
    {
        A_initAsEmpty();
        for (i = begin; i < end; i++) 
        {
            MPI_Recv(&A[i][0], N, MPI_FLOAT, 0, 1, MPI_COMM_WORLD, &status);
        }

    }

    MPI_Barrier(MPI_COMM_WORLD);  //此时每个进程都拿到了数据
    //start_time = MPI_Wtime();
    for (k = 0; k < N; k++) //负责第 k 行的进程进行除法运算，将除法结果一对多广播给编号更大的进程
    {
        if ((begin <= k && k < end)) 
        {
            for (j = k + 1; j < N; j++) 
            {
                A[k][j] = A[k][j] / A[k][k];
            }
            A[k][k] = 1.0;
            for (j = 0; j < total; j++) 
            {
                if(j!=rank)
                    MPI_Send(&A[k][0], N, MPI_FLOAT, j, 0, MPI_COMM_WORLD);//0号消息表示除法完毕
            }
        }
        else
        {
            int src;
            if (k < N / total * total)//在可均分的任务量内
                src = k / (N / total);
            else
                src = total-1;
            MPI_Recv(&A[k][0], N, MPI_FLOAT, src, 0, MPI_COMM_WORLD, &status);
        }
        #pragma omp parallel for num_threads(4)
        for (i = max(begin, k + 1); i < end; i++) 
        {
            for (j = k + 1; j < N; j++) 
            {
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            }
            A[i][k] = 0;
        }
    }
    MPI_Barrier(MPI_COMM_WORLD);	//各进程同步
    if (rank == 0) 
    {//0号进程中存有最终结果
        end_time = MPI_Wtime();
        printf("MPI块划分+openmp耗时:%.4lf ms\n", 1000 * (end_time - start_time));
        //print(A);
    }
}

void xiaoyuan_mpi_openmp_simd(int argc, char* argv[]) //块划分+openmp+simd
{
    double start_time = 0;
    double end_time = 0;
    int total = 0;
    int rank = 0;
    int i = 0;
    int j = 0;
    int k = 0;
    MPI_Status status; //*
    MPI_Comm_size(MPI_COMM_WORLD, &total);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    //cout << rank << " of " << total << " created" << endl;
    int begin = N / total * rank;
    int end = (rank == total - 1) ? N : N / total * (rank + 1);
    //cout << "rank " << rank << " from " << begin << " to " << end << endl;
    if (rank == 0) 
    {  //0号进程初始化矩阵
        A_init();
        //cout << "initialize success:" << endl;
        //print(A);
        //cout << endl << endl;
        start_time = MPI_Wtime();
        for (j = 1; j < total; j++) 
        {
            int b = j * (N / total), e = (j == total - 1) ? N : (j + 1) * (N / total);
            for (i = b; i < e; i++) 
            {
                MPI_Send(&A[i][0], N, MPI_FLOAT, j, 1, MPI_COMM_WORLD);//1是初始矩阵信息，向每个进程发送数据
            }
        }
    }
    else 
    {
        A_initAsEmpty();
        for (i = begin; i < end; i++) 
        {
            MPI_Recv(&A[i][0], N, MPI_FLOAT, 0, 1, MPI_COMM_WORLD, &status);
        }

    }

    MPI_Barrier(MPI_COMM_WORLD);  //此时每个进程都拿到了数据
    //start_time = MPI_Wtime();
    for (k = 0; k < N; k++) //负责第 k 行的进程进行除法运算，将除法结果一对多广播给编号更大的进程
    {
        if ((begin <= k && k < end)) 
        {
            for (j = k + 1; j < N; j++) 
            {
                A[k][j] = A[k][j] / A[k][k];
            }
            A[k][k] = 1.0;
            for (j = 0; j < total; j++) 
            {
                if(j!=rank)
                    MPI_Send(&A[k][0], N, MPI_FLOAT, j, 0, MPI_COMM_WORLD);//0号消息表示除法完毕
            }
        }
        else
        {
            int src;
            if (k < N / total * total)//在可均分的任务量内
                src = k / (N / total);
            else
                src = total-1;
            MPI_Recv(&A[k][0], N, MPI_FLOAT, src, 0, MPI_COMM_WORLD, &status);
        }
        #pragma omp parallel for num_threads(4)
        for (i = max(begin, k + 1); i < end; i++) 
        {
            float32x4_t vaik;
            vaik = vdupq_n_f32(A[i][k]);
            int j;
            for (j = k + 1; j + 4 <= N; j += 4) 
            {
                float32x4_t vakj, vaij, vx;
                vakj = vld1q_f32(&A[k][j]);
                vaij = vld1q_f32(&A[i][j]);
                vx = vmulq_f32(vakj, vaik);
                vaij = vsubq_f32(vaij, vx);
                vst1q_f32(&A[i][j], vaij);
            }
            for (; j < N; j++) 
            {
                A[i][j] = A[i][j] - A[k][j] * A[i][k];
            }
            A[i][k] = 0;
        }
    }
    MPI_Barrier(MPI_COMM_WORLD);	//各进程同步
    if (rank == 0) 
    {//0号进程中存有最终结果
        end_time = MPI_Wtime();
        printf("MPI块划分+openmp+simd耗时:%.4lf ms\n", 1000 * (end_time - start_time));
        //print(A);
    }
}

void xiaoyuan_mpi_circle(int argc, char* argv[]) //循环划分
{  //等步长的循环划分
    double start_time = 0;
    double end_time = 0;
    int total = 0;
    int rank = 0;
    int i = 0;
    int j = 0;
    int k = 0;
    MPI_Status status;
    MPI_Comm_size(MPI_COMM_WORLD, &total);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    //cout << rank << " of " << total << " created" << endl;
    if (rank == 0) 
    {  //0号进程初始化矩阵
        A_init();
        //cout << "initialize success:" << endl;
        //print(A);
        //cout << endl << endl;
        for (j = 1; j < total; j++) 
        {
            for (i = j; i < N; i += total)
            {
                MPI_Send(&A[i][0], N, MPI_FLOAT, j, 1, MPI_COMM_WORLD);//1是初始矩阵信息，向每个进程发送数据
            }
        }

    }
    else 
    {
        A_initAsEmpty();
        for (i = rank; i < N; i+=total)
        {
            MPI_Recv(&A[i][0], N, MPI_FLOAT, 0, 1, MPI_COMM_WORLD, &status);
        }
    }
   
    MPI_Barrier(MPI_COMM_WORLD);
    start_time = MPI_Wtime();
    for (k = 0; k < N; k++) 
    {
        if (k%total==rank)
        {
            for (j = k + 1; j < N; j++) 
            {
                A[k][j] = A[k][j] / A[k][k];
            }
            A[k][k] = 1.0;
            for (j = 0; j < total; j++) 
            { //
                if (j != rank)
                    MPI_Send(&A[k][0], N, MPI_FLOAT, j, 0, MPI_COMM_WORLD);//0号消息表示除法完毕
            }
        }
        else 
        {
            int src = k%total;
            MPI_Recv(&A[k][0], N, MPI_FLOAT, src, 0, MPI_COMM_WORLD, &status);
        }
        int begin = k;
        while (begin % total != rank)
            begin++;
        for (i =begin; i < N; i+=total) 
        {
            for (j = k + 1; j < N; j++) 
            {
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            }
            A[i][k] = 0;
        }
    }
    MPI_Barrier(MPI_COMM_WORLD);	//各进程同步
    if (rank == 0) 
    {//0号进程中存有最终结果
        end_time = MPI_Wtime();
        printf("MPI循环划分耗时：%.4lf ms\n", 1000 * (end_time - start_time));
        //print(A);
    }
}

void xiaoyuan_mpi_circle_simd(int argc, char* argv[])
{  //等步长的循环划分
    double start_time = 0;
    double end_time = 0;
    int total = 0;
    int rank = 0;
    int i = 0;
    int j = 0;
    int k = 0;
    MPI_Status status;
    MPI_Comm_size(MPI_COMM_WORLD, &total);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    //cout << rank << " of " << total << " created" << endl;
    if (rank == 0) 
    {  //0号进程初始化矩阵
        A_init();
        //cout << "initialize success:" << endl;
        //print(A);
        //cout << endl << endl;
        for (j = 1; j < total; j++) 
        {
            for (i = j; i < N; i += total)
            {
                MPI_Send(&A[i][0], N, MPI_FLOAT, j, 1, MPI_COMM_WORLD);//1是初始矩阵信息，向每个进程发送数据
            }
        }

    }
    else 
    {
        A_initAsEmpty();
        for (i = rank; i < N; i+=total)
        {
            MPI_Recv(&A[i][0], N, MPI_FLOAT, 0, 1, MPI_COMM_WORLD, &status);
        }
    }
   
    MPI_Barrier(MPI_COMM_WORLD);
    start_time = MPI_Wtime();
    for (k = 0; k < N; k++) 
    {
        if (k%total==rank)
        {
            for (j = k + 1; j < N; j++) 
            {
                A[k][j] = A[k][j] / A[k][k];
            }
            A[k][k] = 1.0;
            for (j = 0; j < total; j++) 
            { //
                if (j != rank)
                    MPI_Send(&A[k][0], N, MPI_FLOAT, j, 0, MPI_COMM_WORLD);//0号消息表示除法完毕
            }
        }
        else 
        {
            int src = k%total;
            MPI_Recv(&A[k][0], N, MPI_FLOAT, src, 0, MPI_COMM_WORLD, &status);
        }
        int begin = k;
        while (begin % total != rank)
            begin++;
        for (i =begin; i < N; i+=total) 
        {
            float32x4_t vaik;
            vaik = vdupq_n_f32(A[i][k]);
            int j;
            for (j = k + 1; j + 4 <= N; j += 4) 
            {
                float32x4_t vakj, vaij, vx;
                vakj = vld1q_f32(&A[k][j]);
                vaij = vld1q_f32(&A[i][j]);
                vx = vmulq_f32(vakj, vaik);
                vaij = vsubq_f32(vaij, vx);
                vst1q_f32(&A[i][j], vaij);
            }
            for (; j < N; j++) 
            {
                A[i][j] = A[i][j] - A[k][j] * A[i][k];
            }
            A[i][k] = 0;
        }
    }
    MPI_Barrier(MPI_COMM_WORLD);	//各进程同步
    if (rank == 0) 
    {//0号进程中存有最终结果
        end_time = MPI_Wtime();
        printf("MPI循环划分+simd耗时：%.4lf ms\n", 1000 * (end_time - start_time));
        //print(A);
    }
}

void xiaoyuan_mpi_circle_openmp(int argc, char* argv[])
{  //等步长的循环划分
    double start_time = 0;
    double end_time = 0;
    int total = 0;
    int rank = 0;
    int i = 0;
    int j = 0;
    int k = 0;
    MPI_Status status;
    MPI_Comm_size(MPI_COMM_WORLD, &total);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    //cout << rank << " of " << total << " created" << endl;
    if (rank == 0) 
    {  //0号进程初始化矩阵
        A_init();
        //cout << "initialize success:" << endl;
        //print(A);
        //cout << endl << endl;
        for (j = 1; j < total; j++) 
        {
            for (i = j; i < N; i += total)
            {
                MPI_Send(&A[i][0], N, MPI_FLOAT, j, 1, MPI_COMM_WORLD);//1是初始矩阵信息，向每个进程发送数据
            }
        }

    }
    else 
    {
        A_initAsEmpty();
        for (i = rank; i < N; i+=total)
        {
            MPI_Recv(&A[i][0], N, MPI_FLOAT, 0, 1, MPI_COMM_WORLD, &status);
        }
    }
   
    MPI_Barrier(MPI_COMM_WORLD);
    start_time = MPI_Wtime();
    for (k = 0; k < N; k++) 
    {
        if (k%total==rank)
        {
            for (j = k + 1; j < N; j++) 
            {
                A[k][j] = A[k][j] / A[k][k];
            }
            A[k][k] = 1.0;
            for (j = 0; j < total; j++) 
            { //
                if (j != rank)
                    MPI_Send(&A[k][0], N, MPI_FLOAT, j, 0, MPI_COMM_WORLD);//0号消息表示除法完毕
            }
        }
        else 
        {
            int src = k%total;
            MPI_Recv(&A[k][0], N, MPI_FLOAT, src, 0, MPI_COMM_WORLD, &status);
        }
        int begin = k;
        while (begin % total != rank)
            begin++;
        #pragma omp parallel for num_threads(4)
        for (i =begin; i < N; i+=total) 
        {
            for (j = k + 1; j < N; j++) 
            {
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            }
            A[i][k] = 0;
        }
    }
    MPI_Barrier(MPI_COMM_WORLD);	//各进程同步
    if (rank == 0) 
    {//0号进程中存有最终结果
        end_time = MPI_Wtime();
        printf("MPI循环划分+openmp耗时：%.4lf ms\n", 1000 * (end_time - start_time));
        //print(A);
    }
}

void xiaoyuan_mpi_circle_openmp_simd(int argc, char* argv[])
{  //等步长的循环划分
    double start_time = 0;
    double end_time = 0;
    int total = 0;
    int rank = 0;
    int i = 0;
    int j = 0;
    int k = 0;
    MPI_Status status;
    MPI_Comm_size(MPI_COMM_WORLD, &total);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    //cout << rank << " of " << total << " created" << endl;
    if (rank == 0) 
    {  //0号进程初始化矩阵
        A_init();
        //cout << "initialize success:" << endl;
        //print(A);
        //cout << endl << endl;
        for (j = 1; j < total; j++) 
        {
            for (i = j; i < N; i += total)
            {
                MPI_Send(&A[i][0], N, MPI_FLOAT, j, 1, MPI_COMM_WORLD);//1是初始矩阵信息，向每个进程发送数据
            }
        }

    }
    else 
    {
        A_initAsEmpty();
        for (i = rank; i < N; i+=total)
        {
            MPI_Recv(&A[i][0], N, MPI_FLOAT, 0, 1, MPI_COMM_WORLD, &status);
        }
    }
   
    MPI_Barrier(MPI_COMM_WORLD);
    start_time = MPI_Wtime();
    for (k = 0; k < N; k++) 
    {
        if (k%total==rank)
        {
            for (j = k + 1; j < N; j++) 
            {
                A[k][j] = A[k][j] / A[k][k];
            }
            A[k][k] = 1.0;
            for (j = 0; j < total; j++) 
            { //
                if (j != rank)
                    MPI_Send(&A[k][0], N, MPI_FLOAT, j, 0, MPI_COMM_WORLD);//0号消息表示除法完毕
            }
        }
        else 
        {
            int src = k%total;
            MPI_Recv(&A[k][0], N, MPI_FLOAT, src, 0, MPI_COMM_WORLD, &status);
        }
        int begin = k;
        while (begin % total != rank)
            begin++;
        #pragma omp parallel for num_threads(4)
        for (i =begin; i < N; i+=total) 
        {
            float32x4_t vaik;
            vaik = vdupq_n_f32(A[i][k]);
            int j;
            for (j = k + 1; j + 4 <= N; j += 4) 
            {
                float32x4_t vakj, vaij, vx;
                vakj = vld1q_f32(&A[k][j]);
                vaij = vld1q_f32(&A[i][j]);
                vx = vmulq_f32(vakj, vaik);
                vaij = vsubq_f32(vaij, vx);
                vst1q_f32(&A[i][j], vaij);
            }
            for (; j < N; j++) 
            {
                A[i][j] = A[i][j] - A[k][j] * A[i][k];
            }
            A[i][k] = 0;
        }
    }
    MPI_Barrier(MPI_COMM_WORLD);	//各进程同步
    if (rank == 0) 
    {//0号进程中存有最终结果
        end_time = MPI_Wtime();
        printf("MPI循环划分+openmp+simd耗时：%.4lf ms\n", 1000 * (end_time - start_time));
        //print(A);
    }
}

int main(int argc, char* argv[]) 
{
    MPI_Init(&argc, &argv);
    cout <<"问题规模为" <<N <<endl;

    xiaoyuan_cx();
    xiaoyuan_mpi(argc, argv);
    //xiaoyuan_mpi_simd(argc, argv);
    //xiaoyuan_mpi_openmp(argc, argv);
    //xiaoyuan_mpi_openmp_simd(argc, argv);
    //xiaoyuan_mpi_circle(argc, argv);
    //xiaoyuan_mpi_circle_simd(argc, argv);
    //xiaoyuan_mpi_circle_openmp(argc, argv);
    //xiaoyuan_mpi_circle_openmp_simd(argc, argv);
    //mpi_talk(argc, argv);

    MPI_Finalize();
}
